home *** CD-ROM | disk | FTP | other *** search
/ Acorn RISC PD-CD 1 / Acorn RISC PD-CD 1.iso / utilities / _graphics / graphics / _jpeg / c / jrevdct < prev    next >
Encoding:
Text File  |  1991-10-28  |  4.9 KB  |  172 lines
  1. /*
  2.  * jrevdct.c
  3.  *
  4.  * Copyright (C) 1991, Thomas G. Lane.
  5.  * This file is part of the Independent JPEG Group's software.
  6.  * For conditions of distribution and use, see the accompanying README file.
  7.  *
  8.  * This file contains the basic inverse-DCT transformation subroutine.
  9.  *
  10.  * This implementation is based on Appendix A.2 of the book
  11.  * "Discrete Cosine Transform---Algorithms, Advantages, Applications"
  12.  * by K.R. Rao and P. Yip  (Academic Press, Inc, London, 1990).
  13.  * It uses scaled fixed-point arithmetic instead of floating point.
  14.  */
  15.  
  16. #include "jinclude.h"
  17.  
  18.  
  19. /* The poop on this scaling stuff is as follows:
  20.  *
  21.  * Most of the numbers (after multiplication by the constants) are
  22.  * (logically) shifted left by LG2_DCT_SCALE. This is undone by UNFIXH
  23.  * before assignment to the output array. Note that we want an additional
  24.  * division by 2 on the output (required by the equations).
  25.  *
  26.  * If right shifts are unsigned, then there is a potential problem.
  27.  * However, shifting right by 16 and then assigning to a short
  28.  * (assuming short = 16 bits) will keep the sign right!!
  29.  *
  30.  * For other shifts,
  31.  *
  32.  *     ((x + (1 << 30)) >> shft) - (1 << (30 - shft))
  33.  *
  34.  * gives a nice right shift with sign (assuming no overflow). However, all the
  35.  * scaling is such that this isn't a problem. (Is this true?)
  36.  */
  37.  
  38.  
  39. #define ONE 1L                  /* remove L if long > 32 bits */
  40.  
  41. #ifdef RIGHT_SHIFT_IS_UNSIGNED
  42. #define LG2_DCT_SCALE 15
  43. #define RIGHT_SHIFT(_x,_shft)   ((((_x) + (ONE << 30)) >> (_shft)) - (ONE << (30 - (_shft))))
  44. #else
  45. #define LG2_DCT_SCALE 16
  46. #define RIGHT_SHIFT(_x,_shft)   ((_x) >> (_shft))
  47. #endif
  48.  
  49. #define DCT_SCALE (ONE << LG2_DCT_SCALE)
  50.  
  51. #define LG2_OVERSCALE 2
  52. #define OVERSCALE (ONE << LG2_OVERSCALE)
  53.  
  54. #define FIX(x)  ((INT32) ((x) * DCT_SCALE + 0.5))
  55. #define FIXO(x)  ((INT32) ((x) * DCT_SCALE / OVERSCALE + 0.5))
  56. #define UNFIX(x)   RIGHT_SHIFT((x) + (ONE << (LG2_DCT_SCALE-1)), LG2_DCT_SCALE)
  57. #define UNFIXH(x)  RIGHT_SHIFT((x) + (ONE << LG2_DCT_SCALE), LG2_DCT_SCALE+1)
  58. #define UNFIXO(x)  RIGHT_SHIFT((x) + (ONE << (LG2_DCT_SCALE-1-LG2_OVERSCALE)), LG2_DCT_SCALE-LG2_OVERSCALE)
  59. #define OVERSH(x)   ((x) << LG2_OVERSCALE)
  60.  
  61. #define SIN_1_4 FIX(0.7071067811856476)
  62. #define COS_1_4 SIN_1_4
  63.  
  64. #define SIN_1_8 FIX(0.3826834323650898)
  65. #define COS_1_8 FIX(0.9238795325112870)
  66. #define SIN_3_8 COS_1_8
  67. #define COS_3_8 SIN_1_8
  68.  
  69. #define SIN_1_16 FIX(0.1950903220161282)
  70. #define COS_1_16 FIX(0.9807852804032300)
  71. #define SIN_7_16 COS_1_16
  72. #define COS_7_16 SIN_1_16
  73.  
  74. #define SIN_3_16 FIX(0.5555702330196022)
  75. #define COS_3_16 FIX(0.8314696123025450)
  76. #define SIN_5_16 COS_3_16
  77. #define COS_5_16 SIN_3_16
  78.  
  79. #define OSIN_1_4 FIXO(0.707106781185647)
  80. #define OCOS_1_4 OSIN_1_4
  81.  
  82. #define OSIN_1_8 FIXO(0.3826834323650898)
  83. #define OCOS_1_8 FIXO(0.9238795325112870)
  84. #define OSIN_3_8 OCOS_1_8
  85. #define OCOS_3_8 OSIN_1_8
  86.  
  87. #define OSIN_1_16 FIXO(0.1950903220161282)
  88. #define OCOS_1_16 FIXO(0.9807852804032300)
  89. #define OSIN_7_16 OCOS_1_16
  90. #define OCOS_7_16 OSIN_1_16
  91.  
  92. #define OSIN_3_16 FIXO(0.5555702330196022)
  93. #define OCOS_3_16 FIXO(0.8314696123025450)
  94. #define OSIN_5_16 OCOS_3_16
  95. #define OCOS_5_16 OSIN_3_16
  96.  
  97.  
  98. INLINE
  99. LOCAL void
  100. fast_idct_8 (DCTELEM *in, int stride)
  101. {
  102.   /* tmp1x are new values of tmpx -- flashy register colourers
  103.    * should be able to do this lot very well
  104.    */
  105.   INT32 tmp10, tmp11, tmp12, tmp13;
  106.   INT32 tmp20, tmp21, tmp22, tmp23;
  107.   INT32 tmp30, tmp31;
  108.   INT32 tmp40, tmp41, tmp42, tmp43;
  109.   INT32 tmp50, tmp51, tmp52, tmp53;
  110.   INT32 in0, in1, in2, in3, in4, in5, in6, in7;
  111.   
  112.   in0 = in[       0];
  113.   in1 = in[stride  ];
  114.   in2 = in[stride*2];
  115.   in3 = in[stride*3];
  116.   in4 = in[stride*4];
  117.   in5 = in[stride*5];
  118.   in6 = in[stride*6];
  119.   in7 = in[stride*7];
  120.   
  121.   tmp10 = (in0 + in4) * COS_1_4;
  122.   tmp11 = (in0 - in4) * COS_1_4;
  123.   tmp12 = in2 * SIN_1_8 - in6 * COS_1_8;
  124.   tmp13 = in6 * SIN_1_8 + in2 * COS_1_8;
  125.   
  126.   tmp20 = tmp10 + tmp13;
  127.   tmp21 = tmp11 + tmp12;
  128.   tmp22 = tmp11 - tmp12;
  129.   tmp23 = tmp10 - tmp13;
  130.   
  131.   tmp30 = UNFIXO((in3 + in5) * COS_1_4);
  132.   tmp31 = UNFIXO((in3 - in5) * COS_1_4);
  133.   
  134.   tmp40 = OVERSH(in1) + tmp30;
  135.   tmp41 = OVERSH(in7) + tmp31;
  136.   tmp42 = OVERSH(in1) - tmp30;
  137.   tmp43 = OVERSH(in7) - tmp31;
  138.   
  139.   tmp50 = tmp40 * OCOS_1_16 + tmp41 * OSIN_1_16;
  140.   tmp51 = tmp40 * OSIN_1_16 - tmp41 * OCOS_1_16;
  141.   tmp52 = tmp42 * OCOS_5_16 + tmp43 * OSIN_5_16;
  142.   tmp53 = tmp42 * OSIN_5_16 - tmp43 * OCOS_5_16;
  143.   
  144.   in[       0] = UNFIXH(tmp20 + tmp50);
  145.   in[stride  ] = UNFIXH(tmp21 + tmp53);
  146.   in[stride*2] = UNFIXH(tmp22 + tmp52);
  147.   in[stride*3] = UNFIXH(tmp23 + tmp51);
  148.   in[stride*4] = UNFIXH(tmp23 - tmp51);
  149.   in[stride*5] = UNFIXH(tmp22 - tmp52);
  150.   in[stride*6] = UNFIXH(tmp21 - tmp53);
  151.   in[stride*7] = UNFIXH(tmp20 - tmp50);
  152. }
  153.  
  154.  
  155. /*
  156.  * Perform the inverse DCT on one block of coefficients.
  157.  *
  158.  * Note that this code is specialized to the case DCTSIZE = 8.
  159.  */
  160.  
  161. GLOBAL void
  162. j_rev_dct (DCTBLOCK data)
  163. {
  164.   int i;
  165.   
  166.   for (i = 0; i < DCTSIZE; i++)
  167.     fast_idct_8(data+i*DCTSIZE, 1);
  168.   
  169.   for (i = 0; i < DCTSIZE; i++)
  170.     fast_idct_8(data+i, DCTSIZE);
  171. }
  172.